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NOMENCLATURE 


A Coefficient of ())^^ for equation (1) 

B Coefficient of <()^^ for equation (1) 

C Coefficient of (}) for equation (1) 

XX 

C Pressure coefficient 

P 

c Airfoil chord length 

f Slope of the modified airfoil 

f^^ Slope of the original airfoil 

H Shape factor 

I Number of iteration 

i Number of iterations before I 

K Dimensionless frequency 

M Mach number 

m Exponent of coefficient C 

N Transformed coordinate of n 

N^p Number of supersonic points 

n Coordinate normal to the surface 
P Pressure 

R Residual of each iteration 

Re Reynolds number 

S Transformed coordinate of s 

s Coordinate along the surface 

T Temperature 

t Time 

U Transformed velocity of u 

Uoo Free stream velocity 

u Velocity component in the s-direction 
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V Velocity component in the n-direction 

v(s) Surface vertical velocity 

w Wedge thickness 

X Physical coordinate in the flow direction 

y Physical coordinate normal to x 


a 



e 

e 

max 

V 

P 

T 

0 ) 


Angle of attack 

Coefficient for wedge thickness 

Coefficient for vertical surface velocity 

Specific heat ratio 

Airfoil thickness to chord ratio 

Transformed boundary layer thickness 

Displacement thickness 

Momentum thickness 

Maximum wedge angle 

Kinematic viscosity 

Density 

Shear stress 

Velocity potential 

Frequency 


Subscript 

acpt Acceptable value 

e Edge of the boundary layer 

i Ad iaba t ic condition 

max Maximum value 

n Derivative with respect to n 

ref Reference point 
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s Derivative with respect to s 

sh Shock location 

t Derivative with respect to time 

tr Transformed parameter 

w Wall 

X Derivative with respect to x 

y Derivative with respect to y 

0 Stagnation 

1 Upstream of the shock wave 
Free stream 

Superscript 

- Temperature averaged parameter 
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I . INTRODUCTION 


Experimental work by Spaid and Bachalo (1) using holographic inter- 
ferometry produced some clear pictures of the density distribution for 
transonic flow about airfoils. It clearly indicates that there are re- 
gions in which the viscous effect plays an important role: 

1. The shock/boundary-layer intereaction region. 

2. The boundary layer development region. 

3. The wake region. 

In order to predict aerodynamic performances, it is desirable to 
solve the time dependent Navier-Stokes equations as illustrated by Lin 
et al. (2) for low Reynolds number flows. However, the current genera- 
tion of high speed computers has not yet reached the computational speed 
for this approach to be realistic for analyzing transonic flows, as dis- 
cussed by Chapman (3). Steger (4) solved the Reynolds equations for 
transonic flows with a substantial computation time, using a simple clo- 
sure scheme for turbulent stresses. With advancing computer design and 
Improving turbulence modeling, this method has a high potential to be 
the design tool in the near future. However, for practical applications, 
it is still necessary to use viscous correction of an inviscid solution 
for which many efficient computational methods (5,6, and 7) have been de- 
veloped. 

Several approaches are available to consider the viscous effect for 
inviscid analysis. Yoshihara and Zonar (8) used an empirical model of 
viscous ramp to approximate the suddenly thickened boundary layer behind 
a shock wave. It requires the least amount of computation if the contour 
of the ramp can be adjusted so that it will converge to the inviscid so- 
lution. However, the boundary layer effect before the ramp is totally 
ignored. Nash and Scruggs (9) solved the differential boundary layer 
equation for velocity distributions in the boundary layer region. How- 
ever, the sudden pressure increase behind a shock could not be adequately 
treated. Consequently, it becomes very inefficient in determining the 
actual shock location, especially when flow separation appears possible. 
Green, et al. (10) developed a lag-entrainment method to solve the in- 
tegral boundary layer equation together with the lag-entrainment equations 
for the displacement thickness. Collyer and Lock (11) and Melnik, et al. 
(12) applied this method to the inviscid code developed by Jameson (13 and 
14) and obtained some results that agreed with experimental data. However, 
the Green’s lag-entrainment method required a large number of empirical 
constants which have not been physically verified. Moreover, the compu- 
tational time required for the viscous correction made the computational 
efficiency of the inviscid code irrelevant. 

The objective of this study is to develop a viscous correction method 
which not only improves the accuracy of the inviscid solution but also 
maintains its computational efficiency. A viscous ramp (8), can be used 
to partially simulate the suddenly thickened boundary layer behind a shock 
wave. An available method of solving the integral boundary layer equation 
can be used to calculate the displacement thickness. It is noted that a 
conventional integral boundary layer method modifies the sudden increase 
of pressure behind a shock by a lesser pressure gradient with a longer dis- 
tance. Inserting a viscous wedge at the foot of the shock can correct this 
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situation, if the empirical relation for wedge thickness is designed to 
supplement the inadequancy of the boundary layer method. The inviscid 
correction before the shock can be adequately obtained by a conventional 
boundary layer method. The viscous correction behing the shock is ac- 
complished by superpositioning the viscous wedge thickness on the dis- 
placement thickness. For a strong shock situation, the empirical relation 
of the viscous wedge may include the existence of a separation bubble. 

Lee and Van Dalsem (15) developed a viscous correction method for the in- 
viscid full-potential code, TAIR (16) . By comparing with experimental 
data, it gives similar results as other correction methods by improving 
the accuracy of the inviscid solution for moderately strong shock situa- 
tions. However, contrary to other correction methods, it reduces the 
computational time by reducing the numbers of iteration for reaching a 
converged solution. The improved computational efficiency makes this 
method more attractive in correcting inviscid solutions for unsteady state 
maneuvering as well as for three-dimensional wings. 

Inviscid solution for airfoils maneuvering at low frequency unsteady 
motions was obtained by Ballhaus and Goorjian (17) . Rizzetta and Yoshi- 
hara (18) used an order of magnitude analysis to show that the turbulent 
boundary layer of an airfoil reaches steady state during low frequency 
maneuvering. However, the computational time becomes a critical factor 
for practical applications. Owing to the improvement both in accuracy 
and in computational time, this study is to apply the same principle of 
the viscous wedge and the conventional boundary layer for viscous correc- 
tion of an inviscid small disturbances code, LTRAN2 (19). 

II . LTRAN2 


The equation of motion for an unsteady two-dimensional, transonic 
flow with small disturbance assumption, according to Landahl (20; , may 
be written as: 


with 




2B(J) = C(|) + 

xt XX 


yy 


A 

B 

c 




(Y + 1)M„“ (J)^ 


( 1 ) 


Where is the disturbance velocity potential, is the free-stream 
Mach number, y is the ratio of specific heats, and 6 is the airfoil 
thickness-to-chord ratio. The subscripts x,y, and t are the independent 
variables of space and time. The quantities x. v, t and (J) in equation 
(1) have been scaled by c, c/6^/3, amd cd^/^U^o, respectively. For 

an airfoil of chord length c, traveling with a velocity U^q, and executing 
some unsteady oscillatory motion of frequency O), the reduced frequency 
K is defined as: 


K = 0)c/Uoo (2) 

The choice of the exponent m for the coefficient C is somewhat arbirary 
and is taken to be 2 using the Spreiter scaling. At low reduced 
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frequencies, equation (1) can be approximated as; 


which can be obtained 


2B(|) . = C(j) + (}) 
xt ^xx yy 

from the Euler's equation by assuming 
1 - « 1 


( 3 ) 

(A) 


Solutions of equation (3) were obtained using LTRAN2 which was developed 
by Ballhaus and Goorjian (17) for K < 0.2 and modified by Hessenius and 
Goorjian (19) for K < 1.0 with the pressure coefficient defined as: 


Cp = - 2((^^ + K(})^) (5; 

LTRAN2 uses an H-grid which reasonally satisfies the orthogonality condi- 
tion around the airfoil and in the wake region. Using an alternating 
directional implicit (ADI) algorithm and an approximate factorization (AF2) 
scheme, both steady and unsteady solutions can be obtained with reasonable 
computational time for relatively thin airfoils. Consideration of viscous 
corrections needs not only to improve the accuracy but also to maintain 
the computational speed. 

III. SHOCK/BOUNDARY-LAYER INTERACTION 


The occurence of a shock wave not only causes a discontinuity in the 
inviscid flow region but also produces a stronger adverse pressure gradient 
in the boundary layer region. Downstream of a shock wave, the boundary 
layer is suddenly thickened. Sometimes, it may be accompanied by a region 
of flow separation. Theoretically only the Reynolds equations, which are 
obtained from the time averaged Navier-Stokes equations, are adequate in 
analyzing the phenomenon of shock/boundary layer interaction (4). However, 
the required computer time for such a small region makes it impractical for 
industrial applications. Yoshihara and Zonars (8) used a viscous ramp to 
approximate the suddenly thickened boundary layer region behind a shock wave 
with reasonable success. Lee and Van Dalsem (15) developed a method by 
using a simpler viscous wedge superimposed on a turbulent boundary layer 
for correcting steady full-potential solutions. The same principle is being 
applied for unsteady small disturbance solution as follows: 


A. Viscous Wedge 


An empirical formula, which simulates the suddenly thickened boundary 
layer behind the shock, was developed by Lee and Van Dalsem (15) and is 
used here as follows: 


1 0 , 

where 6^ is an empirical constant (B-j = 0.1 


for s < s , 
sh 


s > s , 
— sh 


has been used for both 


( 6 ) 
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full-potential and small disturbance corrections) . Smax is the maximum 
deflection angle for an attached shock at a given upstream Mach number, 
which is determined by the inviscid code, s is the distance along 
the surface of the airfoil with Sg^ as the location of the shock wave. 

For thin airfoils, the suddenly thickened boundary layer can be treated 
as an equivalent vertical surface velocity, v; which is used as a boun- 
dary condition for the inviscid flow solver. 

V (s) = [M (s) • w (s)] (7) 

^ s 

where ^2 empirical constant (p 2 = 2 for the unsteady small distur- 

bance correction) . The vertical surface velocity needs to be scaled by 
6 in order to be treated similarly as other parameters resulted in surface 
variations during unsteady motions. The required additional computing 
time for evaluating the wedge thickness, w, and the vertical surface ve- 
locity, V, is negligible. 

B. Boundary Layer 

Boundary layer always develops before the shock wave. Near the lead- 
ing edge, there is a laminar boundary layer region which is followed by a 
transition region and a turbulent region. In transonic flow, the laminar 
and the transition regions are very small, their effect to aerod 3 mamic per- 
formance is insignificant. A simple method given by Cohen and Roshtko (21) 
was used mainly to provide the initial condition for the turbulent boundary 
layer. Assuming that transition occurs instantaneously at a prescribed lo- 
cation, an integral method for turbulent boundary layer developed by Sasman 
and Cresci (22) was used. The time-averaged continuity and momentum equa- 
tions for two-dimensional, steady, compressible, turbulent flow can be ex- 
pressed as: 


(pu)g + (pv)^ = 0 


( 8 ; 


pu(u) + pv(u) = -(P)^ + (t) 
s n s n 


(9) 


where u and v are the velocity components in the s and n directions, respec- 
tively, and T is the shear stress. Since the temperature variation in the 
boundary layer is significant at transonic speeds, the Mager transformation 
(23) is employed to simplify the integral equation. The transformed coor- 
dinates are: 


■ -/f+) C-)‘ 


(Y+1)/[2(y- 1)] 
ds 


/T \i/2 r 

•it) / t" 


Tg is a function of o; T io the reference temperature. 


( 10 ) 


( 11 ) 
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( 12 ; 


T 

T 

0 


T 1/3 

0.5 w + 0.22 Pr ' + 


(0.5 - 0.22 Pr 


1/3 



o 

where Pr is the Prandtl number. The shape factor, H, and the momentum 
thickness, 6, are related to the transformed shape factor, and the 

transformed momentum thickness, 6^^* as follows: 


pT 1(Y+1)/[2(Y-1)] 

^ ■ kJ 

with H = -g- 


(13) 

(14) 

(15) 


where is the Mach number at the edge of the boundary layer. The dis- 
placement thickness, £*, may then be expressed 


6 * = ( 6 ^^ + 6 *^) 


T \(3y-1)/[2(y-1)J 


ft) 


(16) 


- 


^ Yy+1)/[2(y-D] 


For adabatlc flow, the transformed displacement thickness, 6* , and the trans- 
formed momentum thickness, 6 , are defined as 

’ tr 


6 * 

tr 


■/(-!;)•■ 


^ dN 


(17) 


(18) 


where 6^^ is the boundary-layer thickness in the transformed coordinates, and U 
is the transformed velocity component in the s direction, with 

L/2 


"••ft)' 


(19) 


The transfomed velocity profile In the turbulent boundary layer obeys the 
power law: 


U 

U 



tr 


>i-ll/2 


( 20 ) 
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( 21 ) 


where the adabatlc shape factor is defined as 

6*tr 


- 6 


tr 


The transformed mementum Integral equation can be written as 


de 6 dU 

^ ^ tr e 

dS dS 

e 



(H, ). 

tr 1 



p U 2 

e e 


The transformed moment-of -momentum equation the becomes 


(22) 


d(H ) dU 


1 Ji.i) ^ 

- 11 N \ pV>l 

®tr ' e L 


l(H ) + IjT f /^\ 

PeUe J ytr/ 


(23) 


The shear stress at the wall, is expressed through an empirical rela- 
tion given by Tetervin (24) as 


T 

V 


= 0.123 

-1.561(H ) 


/u 0 ^ r0.268/T Y- 

"• (^) (tJ 


268 


tr 1 


( 24 ) 


where V is the kinematic viscosity at temperature T. The two ordinary 
differential equations. Equations (22) and (23), can be solved simultaneous- 
ly using the 4th order Runge-Kutta method for the shape factor, H, and the 
momentum thickness , 6 . 

A computer code, SLAYER, which was developed by McNally (25), gives the 
solutions of the displacement thickness, 6*, for both the laminar and the 
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turbulent regions. Before the shock wave, the displacement thickness is 
the only source for the vertical component of the surface velocity. After 
the shock wave, the summation of the displacement thickness and the wedge 
thickness becomes the source of the vertical component of the surface ve- 
locity, V. 

V (s) = ^2 [w (s) + 6* (25) 

The integral boundary layer method solves two simultaneous ordinary 
differential equations using relatively small amount of computational time. 
Nevertheless, it is about two orders of magnitude longer than that of the 
viscous wedge. 

C. Wake 

The wake region affects the transonic aerodynamic performance only in 
the vicinity of the trailing edge. Its contribution to drag force can be 
substantial, expecially when separation occurs before the trailing edge. 
Accurate evaluation of the wake region needs to include the possibility of 
flow separation which requires the solution of the Reynolds equations. For 
the interest of industrial application, the constraint of computational time 
calls for a simpler approach. Since the viscous wedge is capable of appro- 
ximating the separation region before the trailing edge, it is advantageous 
to terminate the boundary layer analysis earlier by using a lower value shape 
factor for flow separation. H = 1.8 was used in this study. The wake is 
then considered as an extended airfoil with its effective thickness deter- 
mined as follows: 

i) If no flow separation occurs along the airfoil, the sximmation of the dis- 
placement thickness and the wedge thickness is assumed constant from the 
trailing edge throughout the wake . 

ii) If flow separation occurs along the airfoil, the summation of the dis- 
placement thickness, the wedge thickness and the airfoil thickness at the 
point of separation is assumed constant from the separation point throughout 
the wake. 


IV. METHOD OF SOLUTIONS 

The objective is to consider the viscous effect without substantially 
increasing computational time. It is desirable that the basic algorithm of 
LTRAN2 be modified only when the viscous effect can make a meaningful con- 
tribution. The flow chart in Figure 1 shows both the inviscid and the vis- 
cous solutions for the steady and unsteady calculations. Figure la consists 
of the input information and the grid generation for LTRAN2 . The H-grid, 
once generated, remains unchanged for both steady and unsteady calculations. 
Two possible methods can be used for steady state solutions. The standard 
method, as shown in Figure lb, is an iterative procedure known as AF2 to solve 
the steady state equations. The alternative method, as shown in Figure Ic, 
solves the unsteady equations with the steady state boundary condition. The 
alternative method allows the time interval be chosen as small as necessary 
to obtain a converged solution. Once the steady state solution is known, 
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riGURE lA, FLOW CHART FC3R GRID CONSTRUCTION OF LTRAN2 








FIGURE IB. Flow chart for the standard method to obtain steady state solution 










FIGURE iCe flow chart FQR THE ALTERNATIVE METHOD TO OeTAIN STEADY STATE SOLUTION 






Figure Id shows that the unsteady calculations can be conducted by pre- 
scribing a given function of time for the unsteady motions. 

The viscous correction consists of two major components: the viscous 

wedge, VISCW, and the boundary layer, BLAYER. This section outlines the 
logic and the procedures of the viscous correction. 

A. Logic For Viscous Correction 

The pressure coefficient is closely related to the shock location which 
is sensitive to the viscous wedge approximation. However, the abrupt change 
of the airfoil surface due to the wedge insertion, especially at higher up- 
stream Mach numbers, causes instability for the numerical method. On the 
other hand, the boundary layer displacement thickness before the shock wave 
does not have any significant variations during low frequency unsteady mo- 
tions. Since the computational time required for the wedge thickness is ap- 
proximately 1% of that for the boundary layer development, it is necessary 
to use a minimum number of corrections by the boundary layer method. It ap- 
pears logical to divid the viscous effect into two portions: 

i) The Boundary Layer Portion 

The boundary layer displacement thickness, which is relatively in- 
sensitive to the frequency maneuvering, can be calculated by a conventional 
Integral boundary layer method. It provides the following functions: 

a) The viscous effect, that exists between the leading edge and the shock 
wave, gives a more realistic Mach number upstream of the shock. 

b) A smooth transition for the airfoil contour at the foot of the shock can 
eliminate the need of a precursor for the viscous ramp used by Yoshihara and 
Zonar (8) . 

c) An effective wake thickness can be calculated throughout the wake region. 

ii) The Wedge Portion 

The viscous wedge thickness, which is responsible to the change of shock 
locations at every instant of the low frequency maneuvering, needs to be used 
to complement the boundary layer results as follows: 

a) The suddenly thickened boundary layer behind a shock can be simulated. 

b) The possibility of a separation bubble can be included in the wedge thick- 
ness . 

c) The wake region. can be reasonably approximated even if separation occurs 
before the trailing edge. 

B. Procedures 

Table 1 gives a list of parameters used for the viscous correction of 
LTRAN2 . Both the wedge thickness, w, and the boundary layer displacement 
thickness, 6* are a function of the distance, s, along the airfoil surface 
in the flow direction, x. Different procedures are used for steady and un- 
steady states. 
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PARAMETERS 

CONDITION 

STEADY state. TsQ. 
(IITH ITERATION 

Wt I, SI 

NSPIII :NSP| I -6» V I. I : I. 2 6 

* IS) 

RMAX-iSlO RACPT 

Wl I. SI SWI 1 -1, SI 

RMAX <5 RACPT 

visis-^(m (Wii.si* 

RMAX:^ RACPT 

UNSTEADY STATE 
ITITH TIME STEP 

WIT, SI 


K 1 

V|T.SI»^|m IWIT.SI* i«*(SllJ2 



table I. VISCOUS correction parameters 


i) Steady Calculations 


Assuming the vertical component of the surface velocity is zero, the 
steady calculation for the standard method can be started the same way as 
the inviscid solution. At the (I)th iteration, the wedge thickness, w(I,s), 
can be calculated using Equation (6) when the number of supersonic points, 
NSP(I), remains constant for six consecutive iterations. The boundary condi- 
tion of the inviscid flow solver is then modified by Equation (7) for every 
iteration. The maximum residual, is monitored to compare with an ac- 
ceptable residual, R Is prescribed. The boundary layer displace- 
ment thickness, 6*(st^Sill be evaluated only once when ^ 10 ^acpf 

boundary condition of the inviscid flow solver is then modified oy Equation 
(25) . The values of 6*(s) encounters no significant changes between the lead- 
ing edge and the shock location as approaches to R ^ Since the boun- 

dary layer variation behind the shock is contributed mafn?y by the viscous 
wedge, the wedge thickness, w(I,s), will be recalculated at every iteration 
and is superimposed on the available 6*(s) . When R _< 5 R , the wedge 
thickness is "frozen" so that w(I,s) = w(I-l,s) untif^R ^ . The 

converged steady state solution for viscous correction ?i^then^oBtained . 

Evaluation of this procedure was conducted by comparing with experi- 
mental data as well as with the inviscid solution. The viscous correction, 
generally, gives better agreement with the experimental data and uses less 
numbers of iteration than the inviscid solution for moderately strong shock 
situations. However, there were cases that the standard method failed to 
converge. An alternative method using a smaller time interval allows the 
unsteady procedure to be used for steady calculations. This will be dis- 
cussed in the following section. 

ii) Unsteady Calculations 

At a given time, t, the wedge thickness w(t,s) can be calculated using 
Equation (6). Superimposing on the steady state displacement thickness, 6*(s), 
the boundary condition of the inviscid flow solver can be determined using 
Equation (25) . It is assumed that the low frequency unsteady motion does not 
have any significant effect on the boundary layer development from the leading 
edge to the shock location. The sacrifice in accuracy for the pressure co- 
efficient was found to be insignificant but the savings in computer time was 
substantial. In case that the alternative method was used for steady state 
solutions, the wedge thickness was superimposed on the boundary layer thick- 
ness which was evaluated at a presctibed time step where an approximated 
steady state shock position was established. 


V RESULTS AND DISCUSSIONS 

Both conventional and supercritical airfoils were studied. Thin air- 
foils, such as NACA64A010, can be analyzed using LTRAN2 with reasonable 
accuracy and will not be discussed here. This section considers two super- 
critical airfoils, RAE 2822 and NLR 7301. Their cross-sections are given 
in Figure 2. 
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Figure 2, Investigated Airfoils 
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A. RAE 2822 Airfoil 


The geometry and some experimental data were given by Cook et. al., 

(26) . The RAE 2822 airfoil has a maximum thickness to chord ratio of 
0.1210, The case of Moo * 0.73, a = 3.19, Reoo * 6.5 x 10^ oscillating 1 
degree about x/c = 0.5 at a reduced frequency K • 0.2 was studied. 

i) Steady State Solution 

Due to the large angle of attack, the alternative method of using the un- 
steady state procedure with very small time increment was adopted. Converged 
steady state solutions were obtained at the 2010th iteration, which corresponds 
to a dimensionless time of 17.54, for both the inviscld and the viscous solutions. 
The boundary layer was calculated at the 510th iteration where the shock wave ap- 
proached its final position. Figure 3 shows the Mach number distributions. It 
is noted that the viscous solution gives a slightly lower Mach number before the 
shock. This is due to the boundary layer effect from the leading edge to the 
shock location. The combined effect of the viscous wedge and the boundary layer 
moves the shock closer to the leading edge than the inviscid solution. Figure 
4 shows the pressure coefficient distributions. The viscous correction gives the 
shock location at x/c = 0.49, while the inviscid solution is at x/c = 0.63 in 
comparison with the experimental data at x/c = 0.53. The viscous correction gives 
better agreement with the experimetal data. Nevertheless, it over-corrects the 
shock location. Examining the results of Rizzeta (27) which used the combination 
of the viscous ramp (8) and the Green's lag-entrainment boundary layer method (10;, 
the shock location was at x/c = 0.39, an even higher over-correction than the pre- 
sent method . 

ii) Unsteady State Solution 

Once the steady state solution is reached. The airfoil is forced to execute 
a pitching motion about the midchord, x/c = 0.5, with a reduced frequency K = 0.2 
at an amplitude of = 1.0 degree. The instantaneous angle of attack, a, is 
given as a function of time, t, as follows: 

a = a + a. sin(t) (26) 

o 1 

where a = 3.19 is the steady state angle of attack. Figure 5 shows the compari- 
son between inviscid and viscous solutions at a = 4.18 degrees or t = 39.41. 

Figure 6 shows the same comparison at a = 2.35 degrees or t = 48.12. The viscous 
effect on the magnitude of the pressure coefficient before the shock is relatively 
small. The effect on shock location is significant. To evaluate the unsteady 
results, the magnitude of the pressure coefficient and the phase angle were inte- 
grated between the 4th and 5th cycles. Figures 7 and 8 show the magnitude of the 
integrated pressure coefficient and the phase angle, respectively. The viscous 
effect reduced the maximum magnitude from 45.5 at x/c = 0.67 to 27.8 at x/c = 0.58 
and the phase angle changed from 106/-86 degrees at x/c = 0.77 to 130/-70 degrees 
at x/c = 0.64. The inviscid solution took 40.32 min. CPU time and the viscous 
solution took 40.52 min. CPU time on the AMDHAL 470/v7 computer. It is evident 
that the alternative method for viscous correction improves the accuracy of the 
inviscid solution without any substantial increase of computational time. Moreover, 
in comparison with the more time consuming lag-entrainment method, the conventional 
integral boundary layer method is quite adequate if a simple empirical viscous 
wedge model is being used. 
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UNSTEADY PROCEDURE FOR STEADY SOLUTION 


Figure 3, Mach Number Distribution (RAE 2S22) 
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Figure A, Pressure Distribution (RAE 2822) 



PRESSURE COEFFICIENT 

TRANSONIC SMALL DISTURBANCES ANALYSIS 

WITH VISCOUS CORRECTIONS 

WITHOUT VISCOUS CORRECTIONS 



0*0 0*1 0.a 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 


x/c 

RAE 2822 AIRFOIL MlNF=0.730 ALF=3.19 
UNSTEADY STATE T=»39.41. ALF=4.18DEG. 
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Figure 5, Pressure Coefficient (T = 39.41) 
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Figure 6, Pressure Coefficient (T = 48.12) 
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Figure 7, Magnitude of Pressure Coefficient (RAE 2822) 
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Figure 8, Phase Angle (RAE 2822) 



B. NLR 7301 Airfoil 


The geometry and some experimental data were given by Davis and Mal- 
colm (28) . The NLR 7301 airfoil has a maximum thickness to chord ratio of 
0.1627. Due to the relatively blunt leading edge, Riegel's rule, as sug- 
gested by Rizzeta (27) , was used to modify the slope of the surface contour 
as follows: 

^xm “ ^xo/f^ + <27) 

where f and f are the modified and the original slopes, respectively, 
several cases were analyzed. The case of ^^ 3 o = 0.752, a = 0.37 degree, 

6.21 X 106 oscillating 2.10 degress about x/c = 0.399 at a reduced frequency 
K = 0.4 is discussed in this section. 

i) Steady Solution 

The standard method was used and gave converged steady state solutions. 
The inviscid solution took 467 iterations to converge with 324 grid points 
in the supersonic region. The viscous solution took 130 iterations to con- 
verge with 276 grid points in the supersonic region. Figure 9 shows the Mach 
number distributions. Again, the boundary layer development between the lead- 
ing edge and the shock waves is responsible for the lower Mach numbers before 
the shock, while the wedge and the boundary layer cause the shock position 
to move forward. Figure 10 shows the pressure coefficient distributions in 
comparison with the experimental data. It can be seen that the magnitude of 
the pressure coefficient for the first 60% of the airfoil agrees well between 
the viscous solution and the experimental data. However the shock position 
was under-corrected with an x/c = 0.61 in comparison with the experimental 
data at an x/c = 0.57 and the inviscid solution at an x/c = 0.67. It is noted 
that the viscous wedge and the boundary layer models are exactly the same for 
both the RAE 2822 and the NLR 7301 airfoils. The reason to over-correct one 
and under-correct the other is not readily explainable. 

ii) Unsteady Solution 

Once the solution for steady state is reached, the airfoil is forced to 
pitch 2.01 degrees about x/c = 0.399 with a reduced frequency K = 0.4. The 
instantaneous angle of attack is given as follows: 

a = |sin(t) + k [(^) - ^^^ref 

where “ 0.399, = 2.01 degrees and = 0.37 degrees. The compari- 

son between inviscid and viscous solutions are shown in Figures 11 and 12 for 
a = 1.95 degrees, or t =21.12, and a = 2.43 degrees, or T = 26.35, respective- 
ly, It is noted that when a = 1.95 degrees, the viscous solution gives lower 
pressure coefficient with shock wave moving closer to the leading edge as ex- 
pected, However, when a = 2.43 degrees, the shock location of the viscous so- 
lution is closer to the trailing edge than the inviscid solution. It appears 
that the viscous effect delays the response of the shock position variations 
when the unsteady motion is being executed at a higher frequency. The inte- 
grated results are shown in Figures 13 and 14 for the magnitude of the pressure 
coefficient and the phase angle, respectively. Again to avoid transient ef- 
fect, the integrated results are obtained between the 4th and 5th cycles. It 
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Figure 9, Mach Number Distribution, (NLR 7301) 
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Figure 11, Pressure Coefficient (T = 21.12) 
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Figure 12, Pressure Coefficient (T = 26,35) 
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Figure 13, Magnitude of Pressure Coefficient (NLR 7301) 
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Figure 14, Phase Angle (NLR 7301) 



noted that the viscous solution gives better agreement with experimental 
data than the inviscid solution. Owing to the smaller number of ilterations 
for obtaining the steady state solution, the viscous solution uses less com- 
putational time. The CPU times were 14.85 min. for the viscous solution, 
and 16.87 min. for the inviscid solution on the AMDHAL 470/v7 computer. It 
is evident that the viscous correction not only improves the accuracy of the 
inviscid solution but also reduces its computational time when the standard 
method is used. 


VI CONCLUSIONS 

Modifications of the inviscid small disturbances codes LTRAN2, has been 
completed for considering the viscous effect on airfoils maneuvering at low 
frequency unsteady motions during transonic flight. Assuming the boundary 
layer thickness before the shock does not vary substantially due to low fre- 
quency motions, an empirical model for viscous wedge downstream of the shock 
can be superimposed on a displacement thickness obtained from the conventional 
method of solving integral boundary layer equations . The following conclu- 
sions are reached: 

1. For steady state, the viscous correction requires less computational time 
than the inviscid solution. 

2. For unsteady state, there is no noticeable increase in computational time 
for the viscous correction. 

3. The viscous correction improved the accuracy of the inviscid solution for 
the sutided moderately strong shock situations. 
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